Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :
La bioinfo c'est : <code>MERVEILLEUX</code>.
N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.
Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :
# Je travaille dans mon répertoire de travail qui est associé au projet "dubii2021" sur le cluster de l'IFB.
cd /shared/projects/dubii2021/glelandais/modules-4-5-evaluation/.
# Création de répertoires pour organiser les fichiers.
mkdir raw-data
mkdir quality-controls
mkdir clean-data
Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
cd raw-data
# Comme je vais utiliser les outils en dehors d'un script,
# il me semble que cette commande est importante.
salloc --cpus-per-task=10 --mem=1G
# Chargement du module
module load sra-tools
# Lancement du téléchargement
srun fasterq-dump -p --split-files SRR10390685
# Deux fichiers FASTQ ont été obtenus
Combien de reads sont présents dans les fichiers R1 et R2 ?
# Le nombre de reads est égal au nombre de lignes divisé
# par 4 (un read comporte 4 lignes)
grep -c "^@" SRR10390685_1.fastq
grep -c "^@" SRR10390685_2.fastq
Les fichiers FASTQ contiennent 7066055 reads.
Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
Quelle est la taille de ce génome ?
# Décompression du féichier génome
gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz
# Décompte du nombre de base
cat GCF_000009045.1_ASM904v1_genomic.fna | grep -v "NC_" | tr --delete "\n" | wc -c
# Note : J'ai trouvé de l'aide ici :
# https://dridk.me/genome_chiffre_1.html
La taille de ce génome est de 4215606 paires de bases.
Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
# Téléchargement
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
# Décompression du fichier
gunzip GCF_000009045.1_ASM904v1_genomic.gff.gz
Combien de gènes sont connus pour ce génome ?
# Récupération de la troisième colonne du fichier et
# dénombrement des occurences de gene
cut -f 3 GCF_000009045.1_ASM904v1_genomic.gff | grep -c gene
4536 gènes sont recensés dans le fichier d’annotation.
Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
# Changement de répertoire de travail
cd ../quality-controls
# Chargement du module fastqc
module load fastqc
# Lancement des analyses qualités
srun fastqc ../raw-data/SRR10390685_1.fastq -o .
srun fastqc ../raw-data/SRR10390685_2.fastq -o .
# Rapport MultiQC
module load multiqc
srun multiqc -d . -o .
La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?
car les valeurs de score qualité sont élevées (> Q30) comme le montre le graphique “Per base sequence quality”
Lien vers le rapport MulitQC
Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?
car tous les reads n’ont pas la même longueur
Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
# Utilisation de la formule du cours :
# N*L/G ; avec N = Nombre de lectures, L = Longueur des lectures et G = Taille du génome
# ((7066055 * 2) * 151) / 4215606 = 506.2021
La profondeur de séquençage est de : 500 X.
Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.
# Changement du répertoire de travail
cd ../clean-data
# Chargement du logiciel fastp
module load fastp
# Lancement de fastp
srun fastp --in1 ../raw-data/SRR10390685_1.fastq --in2 ../raw-data/SRR10390685_2.fastq --out1 ./SRR10390685_1_clean.fastq --out2 ./SRR10390685_2_clean.fastq --html ./fastp.html --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json ./fastp.json
Les paramètres suivants ont été choisis :
| Parametre | Valeur | Explication |
|---|---|---|
Ces paramètres ont permis de conserver reads pairés, soit une perte de % des reads bruts.
Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].
Combien de reads ne sont pas mappés ?
reads ne sont pas mappés.
Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
reads chevauchent le gène d’intérêt.
Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.
A work by Migale Bioinformatics Facility
https://migale.inrae.fr
Our two affiliations to cite us:
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France
Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France